#include "Timogram.hpp"

#include <iostream>

#include <windows.h>
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>

#define _USE_MATH_DEFINES
#include <math.h>

#include "MonteCarlo.hpp"
#include "SystemMatrix.hpp"
#include "TimedVoxels.hpp"
#include "Vector3.hpp"
#include "config.hpp"
#include "misc.hpp"

using namespace std;

namespace Dynamic
{

static float image[196 * 128], image1[196 * 128];

ListModeEvent::ListModeEvent() : t_e(0.0), lambda_Le(0.0)
{

}

Timogram::Timogram() : array(nullptr)
{
	this->array = new std::vector<ListModeEvent>[nlors];
}

Timogram::Timogram(const Timogram& v)
{
	this->array = nullptr;
	*this = v;
}

Timogram::~Timogram()
{
	delete [] array;
}

int Timogram::LorIndex(int iphi, int idphi) { return (iphi * (ndets - 1) - iphi * (iphi - 1) / 2 + idphi - 1); }

void Timogram::operator=(const Timogram& m)
{
	if (&m != this) {
		delete[] this->array;
		if (m.array != nullptr) {
			this->array = new std::vector<ListModeEvent>[nlors];
			for (unsigned int i = 0; i < nlors; i++) {
				this->array[i] = m.array[i];
			}
		} else {
			this->array = nullptr;
		}
	}
}

void Timogram::Clear()
{
	for (unsigned int L = 0; L < nlors; ++L) {
		array[L].clear();
	}
}

void Timogram::Zero()
{
	for (unsigned int L = 0; L < nlors; ++L) {
		for (ListModeEvent & lme : this->array[L]) {
			lme.lambda_Le = 0;
		}
	}
}

void Timogram::InitWithOneEvent(double t_e)
{
	this->Clear();

	ListModeEvent lme;
	lme.lambda_Le = 0;
	lme.t_e = t_e;
	for (unsigned int L = 0; L < nlors; ++L) {
		array[L].push_back(lme);
	}
}

void Timogram::randomize(double measurementInterval, int numberOfEvents, SystemMatrix& A_LV, TimedVoxels& x, double activity)
{
	delete[] this->array;
	this->array = new std::vector<ListModeEvent>[nlors];

	double xMaxValue = x.Max();
	int percent = 0;

	cout << "Generating events..." << endl;

	for (unsigned int e = 0; e < numberOfEvents;) {
		// Get a random voxel coordinate for the event
		int vIndex = getRandomVoxel(x);

		// Choose which LOR
		int lIndex = getRandomLorForVoxel(vIndex, A_LV);

		// Woodcock tracking
		double c = activity * 2 / 3;
		try {
			double t_e = woodcockTrackEmissionTime(
			[&c](double t0, double r, double x) -> double {
				double r2 = max(1 - r, 0.001);
				double tRes = t0 - log(r2) / x* c;

				return tRes;
			}, [&activity, &measurementInterval](double t) -> double {
				//return 1.0 / (t + 0.1);

				double sigma = 1.0;
				return Gauss1D(t, sigma) * activity;
				//return t / measurementInterval * activity;
			}, xMaxValue, measurementInterval);

			ListModeEvent lme;
			lme.t_e = t_e;
			this->array[lIndex].push_back(lme);
		} catch (std::exception& ex) {
			// Woodcock tracking was unsuccessful, this emission was out of
			// the measurement interval
		}

		// Move to the next event
		e++;

		int p = e * 100 / numberOfEvents;

		if (p % 10 == 0 && p > percent) {
			percent = p;
			cout << percent << "%" << endl;
		}
	}
}

double Timogram::sumEvents()
{
	double sum = 0.0;
	for (unsigned int l = 0; l < nlors; l++) {
		sum += this->array[l].size();
	}

	return sum;
}

double Timogram::sumLambdaAccumulated()
{
	double s = 0;
	for (unsigned int L = 0; L < nlors; ++L) {
		for (const ListModeEvent & lme : this->array[L]) {
			s += lme.lambda_Le;
		}
	}
	return s;
}

double Timogram::sumLambda(unsigned int lorIndex)
{
	double sum = 0.0;
	for (const ListModeEvent & lme : this->array[lorIndex]) {
		sum += lme.lambda_Le;
	}

	return sum;
}

double Timogram::sumLambda(int iphi, int idphi)
{
	unsigned int lorIndex = LorIndex(iphi, idphi);
	return this->sumLambda(lorIndex);
}

double Timogram::MaxLambdaAccumulated()
{
	double maxlor = 0;
	for (unsigned int L = 0; L < nlors; ++L) {
		double accValue = this->sumLambda(L);
		if (maxlor < accValue) {
			maxlor = accValue;
		}
	}
	return maxlor;
}

double Timogram::MaxLambda()
{
	double maxlor = 0;
	for (unsigned int L = 0; L < nlors; ++L) {
		for (const ListModeEvent & lme : this->array[L]) {
			if (maxlor < lme.lambda_Le) {
				maxlor = lme.lambda_Le;
			}
		}
	}
	return maxlor;
}

void Timogram::PrintSinogram(const char* ofname, double maxvalue, bool bar /*= false*/)
{
	cout << "\n Sinogram: " << ofname << '\n';
	static char buffer[80];
	sprintf(buffer, "%s.tga", ofname);
	FILE* tgafile = fopen(buffer, "wb");
	sprintf(buffer, "%s.ima", ofname);
	FILE* imafile = fopen(buffer, "w");

	int Xres = (bar) ? xres : yres;
	PrintTGAHeader(tgafile, Xres, yres);

	for (int i = 0; i < ndets; i++) {
		for (int j = 0; j < ndets; j++) {
			double alpha = M_PI * (i + 0.5) / ndets;
			double displace = (((double)j + 0.5) / ndets - 0.5) * volume;
			double dphi = 2.0 * acos(fabs(displace) / R);
			double phi = alpha - dphi / 2.0;
			if (displace < 0) { phi += M_PI; }
			double phi1 = phi + dphi;
			if (phi < 0) { phi += 2 * M_PI; }
			if (phi1 >= 2 * M_PI) { phi1 -= 2 * M_PI; }
			if (phi > phi1) {
				double tmp = phi;
				phi = phi1;
				phi1 = tmp;
			}

			int iphi = (int)(ndets * phi / 2.0 / M_PI);
			if (iphi > ndets - 1) { iphi = ndets - 1; }
			int idphi = (int)(ndets * (phi1 - phi) / 2.0 / M_PI + 0.5);
			if (iphi + idphi > ndets - 1) { idphi = ndets + 1 - iphi; }
			double val = this->sumLambda(iphi, idphi) / maxvalue;
			if (idphi <= gap) { val = 0; }
			/*
			double xmin = 0.001;
			double a = -256.0 / log(xmin);

			if (val > 1) val = 1;
			val = -log(val) * a;
			if (val < 0) val = 0;
			if (val > 255) val = 255.5;

			int ival = (int)(val);
			*/

			image[i * ndets + j] = val;
		}
	}

	for (unsigned int i = 0; i < ndets; ++i) {
		for (unsigned int j = 0; j < ndets; ++j) {
			fprintf(imafile, "%f ", image[i * ndets + j] * maxvalue);
		}
		fprintf(imafile, "\n");
	}
	fclose(imafile);

	for (unsigned int i = 0; i < yres; ++i) {
		for (unsigned int j = 0; j < yres; ++j) {
			int I = i * ndets / yres, J = j * ndets / yres;
			image1[i * yres + j] = image[ I * ndets + J ];
		}
	}
	glClear(GL_COLOR_BUFFER_BIT);
	glRasterPos2f(-1, -1);
	glDrawPixels(yres, yres, GL_LUMINANCE, GL_FLOAT, image1);

	sprintf(buffer, "%d", 0); Message(buffer, false);
	if (maxvalue > 1) {
		sprintf(buffer, "%d", (int)maxvalue); Message(buffer, true);
	} else {
		sprintf(buffer, "%4.4f", maxvalue + 0.0005);
		for (int s = 0; s < 5; s++) { buffer[s] = buffer[s + 1]; }
		buffer[4] = '\0';
		Message(buffer, true);
	}
	//      glFlush(); glFinish();
	glRasterPos2f(-1, -1);
	glReadPixels(0, 0, Xres, yres, GL_LUMINANCE, GL_FLOAT, image);


	for (int i = yres - 1; i >= 0; --i) {
		for (unsigned int j = 0; j < Xres; ++j) {
			if (bar && j >= yres - 3 && j < yres) { image[i * Xres + j] = (double)(i - 1) / yres; }
			PrintTGAColor(tgafile, (double)image[i * Xres + j], 0.7);
		}
	}
	fclose(tgafile);
}

void Timogram::InitMonteCarlo(bool Zeromx, MaterialTable& mt, TimedVoxels& x, int Nsamples, double measurementInterval, 
	bool ScatteredOnly /*= false*/, bool Correction /*= true*/)
{
	if (Zeromx) {
		Zero();
	}

	int ScatteredPhotons, TotalPhotons;
	ScatteredPhotons = TotalPhotons = 0;
	double emissionDensity = 10.0;

	for (int i = 0; i < vrows; i++)
		for (int j = 0; j < vrows; j++) {
			//   cout << i*vrows+j+1 << " out of " << vrows*vrows << std::endl;
			Vector3 index(i, j, 0.0);
			double lastEmissionTimePoint = 0.0;

			//int voxsamples = (int) ((double)Nsamples*nvoxels*x.at(i,j)/x.sum + 0.5
			double voxelActivity = x.getAccumulatedValue(i, j);
			double realsamples = (double)Nsamples * nvoxels * voxelActivity / x.Sum();
			int voxsamples = (int)realsamples;
			if (rnd_0_1() < (realsamples - voxsamples)) {
				voxsamples++;
			}
			for (int sampl = 0; sampl < voxsamples; sampl++) {
				Vector3 perturb = index;
				perturb.x += rnd_0_1();
				perturb.y += rnd_0_1();

				int lorIndex = PhotonPairSimulation(index, perturb, mt, voxsamples, ScatteredOnly, Correction, TotalPhotons, ScatteredPhotons);
				
				// If the photon tracked has enough energy
				if (lorIndex != -1) {
					// We start from the voxels last emission time point
					double xvt = x.getValueAtTime(TimedVoxels::VoxelIndex(i, j), lastEmissionTimePoint);
					double T = -log(1 - rnd_0_1()) / (xvt * emissionDensity);
					double emissionTimePoint = lastEmissionTimePoint + T;

					// If the emission time is out of the range of the measurement interval, we are finished with this voxel
					if (emissionTimePoint >= measurementInterval) {
						break;
					}

					ListModeEvent lme;
					lme.t_e = emissionTimePoint;
					lastEmissionTimePoint = emissionTimePoint;
					this->array[lorIndex].push_back(lme);
				}
			}
		}

	cout << "Ratio: " << ScatteredPhotons / (double)TotalPhotons << std::endl;
}

}